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ABSTRACT 

We study mass models that correspond to MOND (triaxial) potentials for which 
the Hamilton-Jacobi equation separates in ellipsoidal coordinates. The problem is first 
discussed in the simpler case of deep-MOND systems, and then generalized to the 
full MOND regime. We prove that the Kuzmin property for Newtonian gravity still 
holds, i.e., that the density distribution of separable potentials is fully determined once 
the density profile along the minor axis is assigned. At variance with the Newtonian 
case, the fact that a positive density along the minor axis leads to a positive density 
everywhere remains unproven. We also prove that (i) all regular separable models 
in MOND have a vanishing density at the origin, so that they would correspond to 
centrally dark-matter dominated systems in Newtonian gravity; (ii) triaxial separable 
potentials regular at large radii and associated with finite total mass leads to density 
distributions that at large radii are not spherical and decline as ln(r)/r 5 ; (hi) when the 
triaxial potentials admit a genuine Frobenius expansion with exponent < e < 1, the 
density distributions become spherical at large radii, with the profile ln(r)/r 3+2e . After 
presenting a suite of positive density distributions associated with MOND separable 
potentials, we also consider the important family of (non-separable) triaxial potentials 
Vi introduced by de Zeeuw & Pfenniger, and we show that, as already known for 
Newtonian gravity, they obey the Kuzmin property also in MOND. The ordinary 
differential equation relating their potential and density along the z-axis is an Abel 
equation of the second kind that, in the oblate case, can be explicitly reduced to 
canonical form. 

Key words: galaxies: kinematics and dynamics — galaxies: structure — dark matter 
— methods: analytical — stellar dynamics 



1 INTRODUCTION 



problems presented by this phenomenological formulation of 
the theory (now known as Modified Newtonian Dynamics 
or MOND), such as conservation of linear momentum (e.g., 
Felten 1984), Bekenstein & Milgrom (1984) substituted the 
heuristic model with the MOND non-relativistic field equa- 
tion 



Milgrom (1983) proposed that the failure of galactic rotation 
curves to decline in Keplerian fashion outside the galaxies' 
luminous body arises not because galaxies are embedded in 
massive dark halos obeying Newtonian gravity, but because 
Newton's law of gravity has to be modified for fields that 
generate accelerations smaller than some characteristic value 
ao ~ 1.2 x 10 -8 cm s -2 . Subsequently, in order to solve basic 




(i) 
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where ||...|| is the standard Euclidean norm, \i is a scalar 
function, <j) is the gravitational potential produced by the 
density distribution p, and 



g = -V0, 



(2) 



is the MOND gravitational field experienced by a test par- 
ticle. For an isolated system of finite mass, eq. (pQ) is sup- 
plemented with the natural boundary condition —> 
for ||x|| — >• oo. Equation (pQ) is obtained from a variational 
principle applied to a Lagrangian with all the required sym- 
metries, so the standard conservation laws are obeyeqj. In 
the regime of intermediate accelerations the function ji is not 
fully constrained by theory or observations, while asymptot- 
ically 

t fort<l, 
for t > 1. 
A common choice is 
t 



(3) 



(4) 



(see also Famaey & Binney 2005, Zhao & Famaey 2006). 

From eq. Q it follows that eq. (pQ) reduces to the Pois- 
son equation when ||V0|| ^> ao, while the limit equation 



V- (||V0||V<£) = 47rGa p, 



(5) 



describes systems for which (or regions of space where) 
||V0|| <C ao, i.e. systems for which the MOND predictions 
differ most from the Newtonian ones. Equation J5}, char- 
acterizing the so-called "deep MOND regime" (hereafter 
dMOND), is then of particular relevance in MOND inves- 
tigations, as this is the regime where one hopes to find the 
most significant differences with the predictions of Newto- 
nian gravity. From the mathematical point of view, the l.h.s. 
of eq. J5} is a special case of the so-called p-Laplace operator 
V(||V0|| p_2 V0). The dMOND case corresponds to p = 3, 
the so-called critical case for 5ft 3 . A large body of mathe- 
matical literature is dedicated to the p-Laplacian, but due 
to its non-linearity it is not surprising that several impor- 
tant questions are still open (e.g., see Lindqvist & Manfredi 
2008). We anticipate that the general results in the present 
paper are independent of the specific form of /x, because they 
just follow from the fact that \i is a function of || V0|| , or they 
are obtained for the dMOND regime. 

A considerable body of observational data seems to sup- 
port MOND well beyond its originally intended field of appli- 
cation (see, e.g., Milgrom 2002, Sanders & McGaugh 2002, 
Famaey h McGaugh 2011), but potential problems of the 
theory have been pointed out by many authors (see, e.g. 
The & White 1988, Buote et al. 2002, Sanders 2003, Ciotti 

1 An alternative non-relativistic formulation of MOND, dubbed 
QMOND, has been proposed (Milgrom 2010), but it is not dis- 
cussed here. 



& Binney 2004, Knebe & Gibson 2004, Zhao et al. 2005, 
Ibata et al. 2011, Galianni et al. 2011). At present, the sit- 
uation is in general considered unsettled. It is thus natural 
to study in detail MOND predictions, in particular focusing 
on dMOND systems, i.e. systems that in Newtonian gravity 
would be dark matter dominated. Unfortunately, MOND in- 
vestigations, especially on the theory side, have been consid- 
erably slowed down by the almost complete lack of aspheri- 
cal density-potential pairs, needed to test the predictions in 
cases more realistic than those described by spherical sym- 
metry. In MOND, the main difficulty to obtain exact as- 
pherical density-potential pairs originates from the fact that 
a simple relation between the Newtonian and the MOND 
gravity fields, produced by an assigned density distribution, 
in general does not exist. In fact, the Newtonian potential 
0n obeys the Poisson equation 



V = 47TG/9, 



(6) 



so that /i(|| V0||/ao)V0 and V0n differ, for assigned p, by 
a solenoidal field S = curlh. In turn, the potential vector 
h depends on p, and so it is apriori unknown. The only 
exception is provided by density distributions in which the 
modulus #n = ||gN|| of the Newtonian field gN = — V0n 
is stratified on surfaces of constant 0n (particularly simple 
cases are those of spherically and cylindrically symmetric 
densities, or densities stratified on homogeneous planes, see 
also Brada & Milgrom 1995; Shan et al. 2008). In such cases 
S vanishes, and 



\\m 

ao 



g gN- 



Equation © coincides with the original MOND formulation 
of Milgrom (1983) and can be solved algebraically for g in 
terms of gN, just by taking its norm. 

A general method to build aspherical and exact MOND 
density potential pairs is presented in Ciotti et al. (2006) 
where, by extending the homeoidal expansion technique in- 
troduced in Ciotti & Bertin (2005) for Newtonian gravity, it 
is shown how a "seed" spherical potential can be deformed 
to an axisymmetric or triaxial shape, leading to analytical 
density-potential pairs satisfying the MOND equation. Un- 
fortunately, the resulting density distributions are not fully 
under control in case of major departures from spherical 
symmetry, and regions of unphysical negative density may 
result. 

For these reasons here we explore a different approach, 
i.e. we focus on the possibility to extend to MOND some 
of the remarkable results obtained in Newtonian gravity for 
potentials separable in ellipsoidal coordinates. This not only 
to better understand the properties of MOND systems, but 
also to develop a new method to generate exact solutions de- 
viating from spherical symmetry for the p-Laplace operator. 
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While we refer to other papers for the full account of separa- 
bility in ellipsoidal coordinates (de Zeeuw 1985b, hereafter 
Z85; de Zeeuw & Lynden-Bell 1985, hereafter ZLB85), here 
we just summarize the properties of Newtonian separable 
potentials in ellipsoidal coordinates relevant to the present 
investigation. Some of these properties are indeed shared by 
similar but non-separable families described in de Zeeuw & 
Pfenniger (1988, hereafter ZP88). Specifically, in Newtonian 
separable systems: 

1. The potential along the long axis of the coordinate 
ellipsoids (the z-axis in the standard convention) is related 
to the density profile along this axis by a linear second order 
ordinary differential equation (ODE). This is the Kuzmin 
property. 

2. For assigned density profile along the z-axis, the ODE 
can be integrated completely. Once the parameters of the 
ellipsoidal coordinate systems are fixed, the solution deter- 
mines uniquely the potential (and so the density) over the 
whole space. The density elsewhere is related to that on the 
z-axis by the so-called Kuzmin formula. Usually, the z-axis 
is the short axis of the density distribution. 

3. The Kuzmin formula shows that the density is every- 
where non-negative if this holds along the short axis. This 
is the Kuzmin theorem. 

4. The Kuzmin formula also shows that density profiles 
that fall off along the z-axis faster than z~ 3 lead to finite 
mass. Density profiles that fall off less steep than z~ 4 become 
spherical at large radii, while for p(z) oc z~ A or steeper the 
models have finite flattening at large radii. In particular, 
density profiles that fall off faster than z~ A lead to density 
distributions that in all other directions falls off as r -4 , so 
that such models are quasi-toroidal (de Zeeuw, Peletier & 
Franx 1986; ZP88). 

Of course, separable potentials are a special - albeit very 
important - class of potentials expressed in ellipsoidal coor- 
dinates. The interest in separable potentials is that, once 
a supporting positive density can be found, then their or- 
bital classification can be done exactly, and equally well 
in MOND or in Newtonian gravity. In addition, as we will 
briefly discuss in the Conclusions, separable potentials may 
allow for a contructive approach towards the assembly of 
self-consistent MOND modes, i.e., collisionless systems sup- 
ported by a positive phase-space distribution function obey- 
ing the Jeans Theorem (e.g., Lynden-Bell 1962, Binney & 
Tremaine 2008). However, when considering the more gen- 
eral problem of obtaining a flexible approach to the con- 
struction of triaxial MOND potential-density pairs, differ- 
ent classes of (non-separable) potentials still obeying the 
Kuzmin property are worth to be explored, such as those de- 
scribed in ZP88. These models are not separable, but their 
simpler algebrical structure leads to simpler equations, that 



might be solved (numerically) more easily than in the sepa- 
rable case. 

The paper is organized as follows. In Section 2 we de- 
rive the ODE along the z-axis for separable models in the 
dMOND regime, and then we extend the result to the full 
MOND case, thus showing that the Kuzmin property holds 
in MOND. In Section 3 we derive the general asymptotic 
behavior at the center and at large radii of the density dis- 
tributions generated by MOND regular separable potentials. 
Some explicit examples of everywhere positive densities as- 
sociated with separable potentials are then presented in Sec- 
tion 4. A discussion of a special family of non-separable 
triaxial potentials (the V\ family of ZP88) is carried out 
in Section 5, where among other findings it is shown that 
the Kuzmin property holds also for V\ systems. The main 
conclusions are summarized in Section 6. In the Appendix 
we list technical details, together with a brief discussion of 
the separable axisymmetric power-law models by Sridhar & 
Touma (1997) in the context of MOND. 



2 THE KUZMIN PROPERTY FOR 
SEPARABLE MOND SYSTEMS 

We consider mass models that correspond to (triaxial) 
MOND potentials for which the Hamilton- Jacobi equation 
separates in ellipsoidal coordinates (A,/i, zy), and investigate 
whether a Kuzmin formula holds for them. This is not ob- 
vious, as the MOND field equation is non-linear and con- 
siderably more complicated than the Poisson equation. We 
recall that the most general form of a separable potential in 
ellipsoidal coordinates can be written as 

, _ FjX) Fiy) ( o, 

9 (A-MA-i/) (p-u)^-X) (v — — fi) 

where F(r) is an arbitrarily function (Z85, ZLB85); the rel- 
evant properties of ellipsoidal coordinates needed in the fol- 
lowing discussion are summarized in Appendix A. 

It turns out that we can obtain information on the 
full MOND case just by restricting to the simpler case of 
dMOND systems, i.e., by focusing on the properties of the 
p-Laplacian. First, we rewrite eq. J5} in the more convenient 
form 

2^l|V0|| 2 



4irGa p= ||V0||V> + 



2||V0|| ' 



(9) 



2 In principle, three different functions Fi(A), i^O-O, and Fs(u) 
may be allowed in separable potentials. Smooth mass models re- 
quire Fi(—a) = F2(—a) and F2(—/3) = F^{—^) together with 
conditions on the derivatives of these functions, so it is in most 
cases no loss of generality to take F\ = F2 = F3 = F (e.g., 
Lynden-Bell 1962, ZLB85). 
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where the explicit expression for the linear differential op- 
erator = (V0, V) in terms of ellipsoidal coordinates is 
given in eq. (|A8[) . The advantage of working with eq. (|9]) 
instead of eq. (|5]) is that one avoids the computation of the 
derivatives of a norm, with the involved square root. 

In principle, by inserting eq. © in eq. ©, and using 
the formulae reported in the Appendix, some heavy algebra 
will give the expression for the density distribution over the 
whole space as a function of F, and its first and second 
order derivative, F' and F" . Not unexpectedly, the resulting 
formula is quite formidable, and of little use. In practice, 
given F(t), the explicit computation can be performed by 
using one of the many available computer algebra packages, 
and some cases will be discussed in Section 4. 

Here instead we are interested in a more general prop- 
erty, i.e. if the Kuzmin formula (or even the Kuzmin theo- 
rem) holds also for MOND separable systems. From eq. (|A2[) 
it follows that the density profile along the whole z-axis of 
a generic density distribution p(A, /x, v) is given by the func- 
tion 

r p(- a ,-/3,r), -l<r< -P 
ip(r)= I p(-a,T,-P), ~P<r<-a (10) 
{p(T,-a,-P), -a<r 

where in the three intervals z 2 = T + 7. In analogy with Z85 
(see also Kuzmin 1956 for the oblate axisymmetric case), we 
now show that the r.h.s. of eq. J9j), evaluated on the z-axis, 
reduces in the three intervals of r to the same second-order 
ODE for the unknown function F(t), thus proving that also 
in dMOND the function ip(r) determines F(t), and so the 
density field everywhere. In other words, the Kuzmin prop- 
erty (and so a Kuzmin-like formula) holds also for separable 
dMOND systems. 

Let us consider the restriction of eq. J9} to the z-axis. 
The Laplace operator satisfies the Kuzmin property, so there 
is nothing to prove, and its expression along the z-axis is 
given in eq. (27) of Z85. An explicit computation then shows 
that the restriction to the z-axis of ||V0|| 2 is given by the 
unique expression 



P>*l|V0|| 2 ], = -4||V0||? [A/> 2(t + i)M\ , 



where 



M 



dN_ _ F"{r) 
dr ~ (r + a)(T + P) 



■ 2F'(r) 



2r + a + 
(r + a) 2 (r + /3) 2 



2F(r) 



3r 2 + 3r(a + P) + a 2 + /T + a/3 



(r + a) 3 (r + /3) 3 
2F(-a) 2F(-P) 



(T + a)*(P-a) (r + P) 3 (P - a)' 



(13) 



+ 



(14) 



Combining the previous results, we obtain the following 
second-order ODE relating, for separable dMOND systems, 
the density profile along the whole z-axis to F(r): 

27rGa ^(r) = VTTi \M\ { [V 2 0], - 2 AT - 4(r + 7)^} -(15) 

Thus, the preliminary result is that the Kuzmin property 
holds in dMOND (i.e., for the p-Laplacian). Furthermore, 
with the aid of the previous results it follows immediately, by 
restriction of eq. (Jj) to the z-axis, that the Kuzmin property 
(and in principle a Kuzmin formula) also holds for the full 
MOND field equation, as the \i function in eq. (Qp depends on 
|| V0||. Unfortunately, the non-linearity of the problem seems 
to prevent the construction of the explicit Kuzmin formula 
even in dMOND, so that the successive analysis of positiv- 
ity of the density as in Newtonian gravity (Z85) cannot be 
performed, and the Kuzmin theorem remains unproven. 

We conclude this general Section by noticing that, as 
pointed out by the Referee, a more general result on Kuzmin 
property can in fact be obtained, encopassing the present 
results and those in Section 5. In practice, with some an- 
alytical work, it can be shown that the Kuzmin property 
certainly hold in Newtonian gravity and in MOND for any 
potential that can be written as a symmetric function of 
ellipsoidal coordinates, i.e., by a function invariant for the 
transformation A — >> \i — » v. 



||V<^=4(t + 7 )A/- 2 , 

in the three intervals spanned by r, where 



M - 



F'(r) 



(2T + a + l3)F( T ) 

(T + a)( T + p) { T + a) 2 { T + py 
F(-a) F(-p) 



+ 



(11) 



(12) 



(r + a) 2 (/?-a) (t + /3) 2 (fi - a) ' 
Remarkably, we note that, with the exception of the factor 



2^/r + 7, no irrationalities are involved in the expression of 
||V0||z. Of course, care is needed in the evaluation of this 
latter quantity, as it contains the absolute value \J\f\ . Finally, 
some algebra shows that the restriction of D^HV^H 2 to the 
z-axis also admits the unique representation 



3 ASYMPTOTIC BEHAVIORS 

In the previous Section we proved that MOND systems 
with separable potentials in ellipsoidal coordinates obey the 
Kuzmin property, and so in principle a Kuzmin formula 
holds for them. Before embarking on the study of the density 
distributions associated with specific separable potentials, 
we focus on the more general question of the asymptotic be- 
havior of the density, both at the center and at large radii 
(for systems with finite total mass). 
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3.1 Behavior at the center 

For a (regular) function F(r) we begin by considering its 
second-order Taylor expansion near the center, i.e. 



F(X) ~ F(-a) + F'(-a)(\ + a) + 
Ffr)~F(-P) + F'(-P)fr + P) + 
F(u) ~ F(- 7 ) + F'(-7)0 + 7) + 



F"( 


-a)(A + a) 2 






F"(- 


-/3)(M + /?) 2 




2 


F"(- 


-7)(^ + 7) 2 


2 



,(16) 



In the expansion above we limit to second order, but the 
present analysis can be carried out (with increasing alge- 
braically complexity) to any desired ordejfl Before embark- 
ing on the following discussion, it is important to recall that 
the function F allows for a linear gauge i.e., two functions 
differing for ar + b, with a and b constants, lead to the same 
function <j> in eq. ©, as can be easily verified by direct sub- 
stitution. With the aid of this gauge it is possible to assign 
two prescribed values to F at two arbitrary points, for exam- 
ple to impose that F(—at) — F(—j) = 0, so that for regular 
F one can write F(r) = (r + a)(r + j)G(t). This form 
has been proved very useful in several investigations (e.g., 
de Zeeuw 1985a, Z85, Hunter & de Zeeuw 1992; Arnold, de 
Zeeuw & Hunter 1994; van de Ven et al. 2003). Here we re- 
frain from using the factorized form, and all the formulae 
are given in full generality: of course, more compact (but 
less symmetric) expressions in terms of the function G can 
be immediately found from those reported here, by simple 
algebraical substitution. 

We focus first on the dMOND regime, recalling that 
near the center A + a ~ x 2 , and analogous relations hold 
for the \i and v coordinates (see eq. [A3]). We consider the 
behavior of the different operators at the r.h.s of eq. (|9]), 
beginning with the Laplacian. As is well known, the applica- 
tion of the Laplace operator to a regular separable potential 
leads to a finite central value 

V 2 0o = (7 - P)F'(-a) + (a - l)F\-p) + a)F\->y) 
2 A 

where 

A = (a-0)(a- 7 )G0- 7 )<O. (18) 

(Z85, eq. [27]). Therefore, barring the special case of null 
derivatives of F at r = —a ,— f3, and —7 (or the more general 
case discussed later), according to eq. @ the value of the 
central density is non-zero in separable Newtonian systems. 

We now move to the term ||V0||, noticing that it ap- 
pears in the denominator of eq. J9}, and so the convergence 
of the density near the origin may be not guaranteed in case 

3 For example, all the results in this Section have been re- 
obtained by using a Taylor series for F(r) truncated at the 10 th 
order (inclusive), and performing the expansions with Mathemat- 
ica. 



of a vanishing gradient left unbalanced by the behavior of 
the term D^||V0|| 2 . Indeed, the vanishing of ||V0|| at the 
center of a regular triaxial potential is expected from geo- 
metrical considerations, and in fact, by using eq. (|16[) we 
find that near the origin 



\\X7<j>\\ 2 ~4(A 2 x 2 + B 2 y 2 + C 2 z 2 ) 



(19) 



where the three constants Ai,Bi, and C\ depend on a, f3 and 
7, and their explicit expression is given in eqs. (|A11[) - (|A12[) . 
Note that the expression above is positive, because the ex- 
panded function is positive definite, and the higher order 
terms cannot affect the sign for sufficiently small displace- 
ments from the origin. In particular, as the three ellipsoidal 
coordinates are independent, the three coefficients must be 
positive, and in fact they are perfect squares. More generally, 
a similar positivity argument holds for the leading term in 
the expansion of ||V0|| 2 independently of the order, i.e., the 
first non-zero term in the expansion is necessarily positive 
near the origin, as we will show in the following. 

We now focus on the last term in eq. J9}. After some 
computation, it is found that near the origin 



£V||V0|| 2 ~ 16CA?x 2 + B\y 2 + C\z 2 ) 



(20) 



where the three coefficients are the same as in eq. (|19[) . 

Finally, by combining the previous results, a simple cal- 
culation shows that in general near the origin 

^ Aj(2Ai + V 2 0o)x 2 + ff 2 (2ffi + V 2 o )j/ 2 + Ci(2d + V 2 o )z 2 
9 ~ 2nGa ^ A\x 2 + B\y 2 + C 2 z 2 

A change to spherical coordinates then proves that the cen- 
tral density of dMOND separable systems vanishes, linearly 
with the spherical radius r. 

A natural question arises, i.e. can we say something 
about p in the special circumstance of Ai — Bi = C\ — 0? 
It could be that the order balance between the numerator 
(17) and the denominator in eq. J9} breaks down when the lead- 
ing term of ||V(/>|| 2 near the origin is of higher order. The 
obtained results are indeed quite interesting. First of all, by 
inspection of eqs. (|A11[) - (|A12[) . it follows that the condition 
A\ — B\ — Ci = is equivalent to the requirement that 
F f (— a), F f (—/3), and F'(— 7) are well defined functions of 
a, f3, 7, and of F(—a), F(—f3), F(— 7). Accordingly, we con- 
sider the potential in eq. (|16|) . with the values of F' fixed by 
the special case just described (and increasing the adopted 
order of expansion of F to the third one) . The leading terms 
of the expansions now read 

V 2 - 6(^ 2 x 2 + B 2 y 2 + C 2 z 2 ), 

||V0|| 2 ~ 4(A 2 2 x 6 + B 2 y 6 +C 2 z 6 ), (22) 
P ||V0|| 2 ~ 48(^|x 8 + B 3 2 y 8 + C 2 3 z 8 ), 

where the explicit form of the coefficients is given in 



eqs. (|A13|) - (|A14[) . Note that in this case also the Laplace op- 
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erator vanishes at the orgin, and a simple calculation shows 
that the density near the center vanishes as r 5 . 

We finally repeat the argument above, with the addi- 
tional request that also Ai — B2 = C2 = 0, thus fixing the 
values of F"(-a), F"(-(3), and F"(- 7 ) from eqs. (lAl3l) - 
()A14|) . By using eq. (|16|) with F(r) expanded up to the 
fourth order inclusive, we now find 

V 2 <j>~5(A 3 x A + B?,y A +C 3 z A ), 

\\V4>\\ 2 ~Ax w + Biy™+Clz w , (23) 
X^||V0|| 2 ~ lO^lx 14 + Bly 14 +C3V 4 ), 

where the coefficients are given in eqs. (|A15[) : now the den- 
sity at the center vanishes as r 9 . By repeating this explo- 
ration to higher and higher orders, we find that the order 
of vanishing of the density, as a function of radius r, in- 
creases by 4 for each additional order of regularity imposed 
on || V0|| near the origin. We also found that, from the third 
order upward, the first non-zero coefficients Ak, Bk, and Ck 
of the leading terms near the origin depend only on the 
derivatives F^ k \r) evaluated at —a, —f3 and —7, so that 
high-order regularity at the center can be expressed directly 
as the vanishing of the corresponding derivatives of F at the 
origin. 

Therefore, the previous analysis shows that a generic 
dMOND system associated with a regular triaxial separable 
potential has - at variance with Newtonian gravity - a zero 
density at the center. This fact leads to some non-trivial 
consequences. The first is that this result remains true even 
when using the full MOND equation, as can be verified with 
a formal expansion. A simple argument is as follows. If we 
use the full MOND equation, and the central regions are not 
in dMOND regime, then they are described by Newtonian 
gravity. But the Newtonian force in triaxial regular separa- 
ble potentials vanishes, so all MOND separable systems at 
the center are actually in dMOND regime, with the conse- 
quent vanishing central density. Of course, the vanishing of 
the central density in MOND may well occur also for other 
families of non-separable potentials (e.g., for potentials with 
a sufficient degree of reflection symmetries along the coor- 
dinate axes). 

The second consequence follows from a further argu- 
ment. Suppose MOND holds, and consider a separable sys- 
tem of baryonic density p, so that from the previous result 
p = at the origin. We now focus on the the total den- 
sity /On of the so-called Equivalent Newtonian System associ- 
ated with the baryonic density p, i.e., the mass distribution 
needed in Newtonian gravity to produce the same gravity 
field of MOND. Of course, />n is obtained by application of 
the Laplace operator to the MOND potential, so that in the 
Newtonian framework the baryonic density p results "im- 
mersed" in a dark matter halo of density ph = pN — p, and 
from eq. (17) the halo density at the center will be different 



from zero. Provided ph > everywhere (a non trivial re- 
quest), we are lead to conclude that a separable MOND sys- 
tem would appear, when interpreted in the context of New- 
tonian gravity, fully dark matter dominated at the center, 
with an arbitrarily large (formally infinite) mass-to-light ra- 
tio near the origin. Note that this property may be expected 
also in other families of MOND potentials, not necessarily 
separable. In fact, it can be easily proved that a regular 
potential with reflection symmetries - <j) = 0(x 2 , y 2 , z 2 ) - 
leads to a density with an expansion near the center iden- 
tical to eq. (21), independently of separability. However, at 
higher orders the special form of eqs. (22)- (23) is not ob- 
tained, showing that separability removes the cross-terms 
and leaves diagonalized quart ics, sextics, and so on. 

We conclude by noting that the addition of a central 
black hole would change the central force field from dMOND 
to Newtonian, in principle opening the possibility to have 
systems with a non-zero central density. Unfortunately, the 
addition of a central mass breaks down separability of tri- 
axial potentials (excluding exceptional axisymmetric cases, 
such that discussed in Appendix B). 



3.2 Behavior at large radii 

The other place where the asymptotic analysis can be carried 
out in generality is at infinity. In particular, from eqs. (|A1[) - 
(|A2[) , it follows that A ~ r 2 for r — >• 00, with r being the 
radius in spherical coordinates. In order to better illustrate 
the MOND case, we begin by recalling the idea behind the 
computation in Newtonian gravity. For a system of finite 
total mass M, Newtonian gravity dictates that (f) ~ —GM/r 
for r 00. If we ask also for separability, then it is easy 
to show that the required asymptotic trend is matched in 
eq. JHJ) if and only if 

F[t) = r 3/2 /i(r), /i(r)~l for r -> 00. (24) 

Note that in the expression above we fixed GM — 1: it is 
simple to restore this dimensional factor in the obtained den- 
sity distribution after the application of the Laplace opera- 
tor. From now on we refer to h(r) as to the shape function: 
its relevance in determining the mass profile at large radii 
will be discussed in detail in the dMOND context. 

By evaluating the Laplace operator, one recovers the 
well known result of Newtonian gravity (e.g. de Zeeuw, 
Franx & Peletier 1986) that regular separable potentials in 
ellipsoidal coordinates, associated with finite total mass and 
finite flattening at large radii, lead to density distributions 
that share the asymptotic radial behavior (in general mod- 
ulated by angular dependence) 

p oc -i-, r -> 00, (25) 
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with additional properties listed in Point 4 in the Introduc- 
tion. 

We now use a similar approach in MOND. For a fi- 
nite mass system, the leading monopole term of the MOND 
potential follows from eq. jTj), with ~ \/GMao lnr, due 
to the r _1 decline of the MOND acceleration at large dis- 
tances (where the field is weak, and so the system is actually 
dMOND). Therefore, if we ask for separability we are now 
forced to assume 

r 2 h(r) lnr 



F(r) 



h(r) 



1 for r 



(26) 



where the coefficient 1/2 takes into account the relation A ~ 
r 2 . Again, the dimensional coefficient y/GMao is set equal 
to 1: in the density profile obtained by the application of 
the MOND operator, due to its non-linearity, the resulting 
coefficient is GMao. 

As in the Newtonian case, the detailed behavior at in- 
finity of the shape function h(r) determines the radial profile 
of the density distribution (note that this is not quite the 
same as the three-dimensional shape, which is also a func- 
tion of the angular coordinates). In fact, also the first and 
second derivatives of h are involved in the computation of 
p and, as is well known, asymptotic properties of functions 
in general are not shared by their derivative^. In particu- 
lar, an arbitrary choice of h can lead to a system with neg- 
ative densities or infinite total mass, in contrast with the 
hypothesis behind eq. (|26j) . For this reason, the convergence 
of the total mass must be checked for any specific choice 
of h(r). In order to carry out a sufficiently general analysis 
(encompassing a large fraction of the cases arising in prac- 
tical situations), here we restrict to shape functions h with 
a Frobenius expansion for r — » oo, i.e. we assume 
B 



h{r) 



T £ T 1 



+ o 



( T 2+ £ ) 



for t 



(27) 



with e > 0. Note that eq. (|26|) imposes A = for e = 0, while 
for e integer we are actually dealing with a regular function 
at infinity. 

We begin with the regular case e = 1. The computation 
of the dMOND operator does not pose special difficulties in 
the regular case, and for A — ► oo one finds 

1 2(a + /3 + 7 + /x + z/-A/2)lnA 
A 



l|V0|| 



A 2 

2 1 2(A + /i + z/)lnA 



X 



A 2 



(28) 



^l|V0|| 2 



2 | 10(A + /x + i/)lnA | _ 3) 



A 2 A 3 
so that, after restoring the dimensional factor, we obtain for 

A —> oo 

4 An elementary example of a function asymptotic to 1 with ar- 
bitrarily large derivative for x — > oo is 1 + x~ x sinx 3 . For a dis- 
cussion of the possible issues involved in the differentiation of 
asymptotic relations, see e.g. Bender &; Orszag (1978). 



p ~ M (p + u + 4A-2a-2P-2,)lnX + Q (a _ 5/2) ^ 

A few important points should be noted. First, only the 
leading coefficient A appears, while all the higher order co- 
efficients do not affect the leading term of the density ex- 
pansion: at infinity, systems with A = are indistinguish- 
able from systems with constant h = 1. Second, the density 
distribution at large radii is not spherically symmetric, a 
consequence of the exact cancellation of the leading terms 
in eq. (|28p when combined in eq. (0. Third, at large radii 
the density is nowhere negative for A > (2a + 3/3 + 3^)/ 4: 
this posit ivity result is not expected a priori, especially when 
considering the non-linear nature of the p-Laplacian. Finally, 
the radial behavior of the density at large radii is propor- 
tional to ln(r)/r 5 . Curiously, the light distribution of ellip- 
tical galaxies in their external regions seems to be described 
better by the 1/r 4 profile (e.g., see Jaffe 1983, Bertin & Sti- 
avelli 1984, 1989; Hernquist 1990, Dehnen 1993, Tremaine 
et al. 1994, see also Bertin & Stiavelli 1989), characteristic 
of the regular Newtonian separable case. 

The discussion above concludes the case of an h function 
with a regular expansion at infinity. Moving to the case of 
non-integer e, with some additional work it can be shown 
that for e > 1 eq. (|29|) still holds with A = 0, thus extending 
to the irregular cases the result on the effect of higher order 
terms obtained for a regular function h. Therefore, we are 
left with the irregular case < e < 1, i.e. when the shape 
function h(r) admits a genuine Frobenius expansion at oo. 

Lengthy algebra and a careful order balance show that 
the following rigorous formulae, extending the validity of 
those in eq. (|28|) , hold: 



1 _C_ 

A + A!+ e 



2(a + P + 7 + M + v) In A + ^-2 ) 



l|V0|| 2 



A 2 

(1 + D/A ) 2 2( M + i/)lnA + 2 



^l|V0|| 2 



A A 2 
2(l + L>/A e ) 2 (l + £/A' 



(30) 



+ 10( M + z,)lnA +(9(A _3 ) 



A 2 A 3 
where 

C = A[l - 4e - e(l - 2e) In A], 

D = A(l-elnA), (31) 
£ = A[l+4e-e(l + 2e) In A]. 

For e = 1 we recover the results in eq. (|28[) . Combining the 
expansions above in eq. J9]), and expanding the norm at the 
denominator, some additional work finally shows that 

M4Ae(elnA-2)(l + C/A 



4tt A 3 / 2 + e + 

M (p + v - 2a - 2f5 - 2 7 ) In A 



-5/2\ 



(32) 



4tt A 5 / 2 
and again for e = 1 the regular case in eq. ([29]) is recovered. 
The density profile at large radii consists of the Frobenius 
contribution dependent on A, e, and A, and in the regular 
part, independent of e. The first component becomes spher- 
ical at large radii, at variance with the regular component, 
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Figure 1. Left panel: density profiles (normalized to M/4tt, where M is the total mass of the model) along the Cartesian coordinate 
axes, at radial distance r from the center (x: solid, y: dotted, and z: dashed) for a separable model with F(r) = — r 2 ln(r)/2 and a = —3, 
P = —2, and 7 = —1. Note how the density vanishes at the center. Right panel: Density ratios p(x, 0, 0)/p(0, 0, z), p(0, y, 0)/p(0, 0, z), and 
p(x, 0, 0)/p(0, y, 0) as a function of distance from the origin, i.e. x = y = z = r, where p is intended expressed in Cartesian coordinates. 
Lines are, in order, solid, dotted, and dashed, respectively. 




x x y 



Figure 2. Isodensity contours in the three Cartesian coordinate planes for the reference model, given by F(t) = -t 2 ln(r)/2 (see Fig. 1). 
Darker gray corresponds to lower values of the density. From the figure it is apparent how the density decreases near the center and at 
large radii, so that the model presents an ellipsoidal corona of higher density. Note also how the z axis corresponds to the shortest axis 
of the density distribution in the radial interval shown. 



dependent also on the jjl and v coordinates. In addition, 
the spherical component is dominant for < e < 1, being 
asymptotic to AAe 2 ln(A)/A 3/2+e . Therefore, when < e < 1 
the density at large radii is spherical and positive for A > 0, 
while for e > 1 it is non-spherical: it is always positive for 
e > 1, and for e = 1 positivity is assured for A greater than 
some negative value. We conclude by noticing that the choice 
of the potential ([26]) implicitly assumes a finite total mass, 



and in fact this is found in the solution, as p oc ln(r)/r 3+2e 
for < e < 1 and p oc ln(r)/r 5 for 1 < e. 

4 EXPLICIT CASES 

From the general analysis in Section 3 we found that the 
central density of MOND models vanishes for regular po- 
tentials separable in ellipsoidal coordinates. We also found 
that at large radii the positivity is assured, provided a cer- 
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tain coefficient in the series expansion of the shape function 
at infinity is larger than a threshold value (0 in a genuine 
Frobenius expansion, and negative in the regular case, with 
the specific value given by a linear combination of the three 
axial coefficients a, f3 and 7 defining the ellipsoidal coordi- 
nate system). Therefore, we have at least indications that 
everywhere positive triaxial densities with separable poten- 
tials may exist in MOND. 

Unfortunately, without the explicit Kuzmin formula, in 
general the positivity of the density associated with an as- 
signed F(t) can be checked only numerically. We now show 
that such cases in fact can be found routinely. For illustra- 
tive purposes, we start with a representative dMOND posi- 
tive separable model, then we illustrate with a few examples 
how the shape function h affects the resulting densities. For 
simplicity, all the examples presented in this Section are con- 
sidered in the dMOND regime. The use of the full MOND 
equation would introduce a dependence on total mass, lead- 
ing to a more complicated discussion, without adding new 
information to the present discussion. 

The reference model is perhaps the simplest possible, 
and it is obtained by using F(r) — — r 2 ln(r)/2 in eq. (j8|). 
i.e., we just fix h = 1 in eq. {26]) . We also assume a = —3, 
P = — 2, and 7 = — 1, but we stress that global positivity has 
been found for all the explored values of the three parame- 
ters, even in the cases characterized by very different values 
of the axial parameters. The relations between ellipsoidal 
and Cartesian coordinates needed for the plots of the iso- 
density contours in the three coordinate planes are reported 
in Appendix A3. 

In the left panel of Fig. 1 we show the density pro- 
files p(x, 0, 0), p(0, y, 0), and p(0, 0, z) of the reference model, 
where it is intended that now the density is expressed in 
Cartesian coordinates, and x — y = z — r, where r is 
the radial distance from the origin. It is apparent how the 
density vanishes at the center, then reaches a peak, and fi- 
nally declines again, in accordance with the previous asymp- 
totic analysis. The model is not spherical, as can be seen 
from Fig. 2, where we present the density cross-sections 
(in the central regions) in the three Cartesian coordinate 
planes: dark grays correspond to low density values, while 
light grays are the density peaks. All the main features of 
Fig. 1 can be easily recognized, and in particular the density 
depression in the inner regions: overall the density of the 
reference model is characterized by a very nice ellipsoidal 
shape. In practice, the resulting density distribution looks 
similar to an heterogeneous ellipsoid with a non-monotonic 
density stratification. The radial trend of the density shape 
is quantified, on a much larger radial interval, in the right 
panel of Fig. 1, where the density ratios p(x, 0, 0)/p(0, 0, z), 
p(0, y, 0)/p(0, 0, z), and p(x, 0, 0)/p(0, y, 0) are represented 



by the solid, dotted, and dashed lines, respectively, for 
x — y — z — r. Note that in general the density ratios 
p(r, 0, 0) / p(0, 0, r) , and p(0, r, 0) / p(0, 0, r) , for a density dis- 
tribution flattened along the z direction, are > 1 (< 1) if 
the density is decreasing (increasing) with r. A visual in- 
spection of Fig. 2 then explains the behavior of the density 
ratios up to r ~ 100: in these regions the z axis (the long 
axis of the coordinate ellipsoids) corresponds to the "short" 
axis of the density distribution, thus confirming that this 
property usually holds not only in Newtonian separable sys- 
tems, but also in MOND. We also note the interesting oc- 
currence of a double switch between the intermediate and 
the long axis. However, for large r, the density ratios drop 
again below 1, i.e., the density distribution in these regions 
(not shown in Fig. 2) becomes elongated along the z axis: it 
has been verified that this non-sphericity is in perfect agree- 
ment with eq. (|29j) evaluated for A = 0. Therefore, this 
model represents a counterexample (in MOND) for the Ed- 
dington (1915) conjecture, fully discussed in by de Zeeuw et 
al. (1986, Section 4 therein). 

We stress that the remarkable ellipsoid-like density dis- 
tribution of the reference model (preserved also for signifi- 
cantly different values of the axial parameters a, j3, and 7) is 
not a general property of MOND models with separable po- 
tentials. In models where we allow for a non-constant shape 
function, the resulting (positive) densities are quite peculiar, 
in some cases with high-density, detached lobes along the z 
axis or, in other cases, by the presence of curious low-density 
regions. Of course, in accordance with the asymptotic anal- 
ysis, the central density of all these models still vanishes. In 
Fig. 3 we show a suite of such densities obtained for different 
choices of the function h, listed in the caption, and for the 
same values a — — 3, /3 = — 2, and 7 = —1 of the reference 
model in Fig. 2. In particular, moving from the top to the 
bottom rows the coefficient A in eq. (27) decreases, and the 
corresponding densities become more and more complicated. 
This is not surprising, because the coefficient A approaches 
the positivity limit for the density at large radii discussed af- 
ter eq. (|29p . Finally, the comparison of the model in the last 
row with the reference model in Fig. 2 shows how higher- 
order terms in the expansion of h may affect the density, as 
A = in both cases. 



5 THE Vi FAMILY OF DE ZEEUW & 
PFENNIGER (1988) 

As discussed in the Introduction, this paper focuses on 
MOND triaxial models with the potential separable in el- 
lipsoidal coordinates. Separable models are a very special 
subset of triaxial models, and so it is of some interest to see 
what properties of separable models are in fact shared by 




Figure 3. Isodensity contours in the coordinate planes for different choices of the shape function h in eq. (|26|) , and for a = —3, (3 = —2, 
and 7 = -1. From top to bottom: h(r) = 1 + 4/r + 11/r 2 , 1 + 3/r, 1 + 2/r + l/2r 5 , and 1 + 1/r 3 - \/2/r 6 . Dark grays correspond to 
lower density values. 



more general triaxial models, and what are the main prop- 
erties of the density distributions associated to some MOND 
(non separable) potentials in ellipsoidal coordinates. Here we 
restrict for sake of simplicity to the important Vi family 



F(X) + F(n) + F(h 



(33) 



introduced and fully discussed in ZP88. V± potentials are 
relevant here because in Newtonian gravity they obey the 
Kuzmin theorem; in addition, altough non-separable, they 
are algebraically simpler than separable potentials (which 
comprise the family V2 and can be derived from Vi by ap- 
plication of a linear operator [ZP88]). Based on the results 
obtained for separable models, it is reasonable to expect that 
also Vi potentials in MOND satisfy the Kuzmin property. 

In fact, we now show that the restriction to the z-axis 
of the r.h.s. of eq. (9} with the potential (|33|) . reduces in 
the three intervals of r in eq. (|10f) to the same second-order 
ODE. This proves that the Kuzmin property holds also for 
Vi potentials in dMOND regime (i.e., for the p-Laplacian). 
The Laplace operator satisfies the Kuzmin property, so there 
is nothing to prove, and its expression along the z-axis is 
given in Sect. 3.1 of ZP88. An explicit computation then 
shows that over the whole z-axis 



||V0||^=4(r + 7 )F'(r) 



and 



[^IIV^H 2 ], = 4||V0||* [F'(t) + 2(t + 7)f"(r)] 



(34) 



(35) 



The following second-order ODE along the z-axis for 
dMOND Vi systems is finally obtained: 



7r Gao^(r) = (T + 7) 3/2 |F / (r)|x 
4F"(t) + 
(7 



T + 7 

a)F'(-a) 



+ 



+ ■ 



r + a ' r + P 
( 7 -0)F'H3) 



F'(r) 

(r + 7)(r + a) (r + 7 )(r + p) J ' ^ 

Again, in analogy with the separable case, it is not diffi- 
cult to show that the Kuzmin property (and in principle 
a Kuzmin formula) also holds for Vi potentials in the full 
MOND regime. 

As expected, eq. (|36|) is considerably simpler than the 
corresponding eq. JT5J), and some additional classification 
and elaboration can be carried out. In fact, albeit eq. (|36|) 
is still second order, non-homogeneous and non-linear, the 
function F(r) is missing, so that for general Vi systems in 
dMOND, the z-axis ODE can be reduced to a non-linear first 
order equation, solving for F' . In particular, in each r inter- 
val where the sign of F' is constant, the ODE belongs to the 
important family of Abel differential equations of the second 
kind: 

[go (x) + gi (x)y]y' = f (x) + fi(x)y + f 2 (x)y 2 + f 3 (x)y 3 , (37) 
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(e.g., Kamke 1948, Zwillinger 19970, a generalization of the 
Riccati equations (Ince 1964). Clearly, the problem is com- 
plicated by the fact that the sign of F' is not known a pri- 
ori: in the following discussion we assume, for simplicity, 
that F' does not change sign for r > —7, and so we set 
F'(t) — =biJ(r), with H(r) > 0. Under this assumption, 
eq. (|36p can be rewritten as 

HH' = ±f ( T ) + f 1 (T)H + f 2 (T)H 2 , (38) 

so that in our case we have an Abel equation with go = fa — 
0, gi = 1, and 



/o 

h 
h 



nGao^ir) 
4(r + 7) 3 / 2 ' 

( 7 -a)H(-a) , (7 -/?)#(-/?) 
4(r + 7)(r + a 

1 ( 2 
~4 



+ 



+ 



4(r + 7 )(r + /3)' 
1 1 



(39) 



+ 



r + 7 ' r + a ' r + /3 y 

The choice of the sign in front of fo determines (if they 
exist) two solutions H±(r) > 0, so that the problem is fi- 
nally solved for F'(r) — ±H±(r). Unfortunately, the gen- 
eral solution of Abel ODEs (first and second kind) is not 
known, but several remarkable transformations have been 
found (e.g., Polyanin & Zaitsev 2003). For example, by set- 
ting H = g(r)p(r) it is possible to determine p(r) so that 
eq. ([38)) can be written in the reduced (but not yet canoni- 
cal) form 



gg' = Rg±S, R= f -, S=^r 



fo 



V 



m our case 



p(r) = 



(40) 



(41) 



y/TTl\r + a| 1 /4| r + / 3|i/4- 

The canonical form would be then obtained by requiring 
R{t) — 1, through the definition of the new independent 
variable £ = j Rdr. Unfortunately, in the triaxial case the 
evaluation of the integral requires Appell functions, so that 
the inversion r — t(£) is impossible in closed form. In the 
prolate case (ft = 7) the integral reduces to a standard hy- 
pergeometric function, and inversion is again impossible, but 
in the oblate case (ft = a) 



£=(a- 7 )ff(-a)< 



7 + T 



7- 

arcsinh 



-7 < T < 



a + r 



(42) 



7" 



-ol < r: 



and so the reduction to the canonical form is possible by us- 
ing circular and hyperbolic functions. However, the resulting 
Abel equation is still unsolvable. 

As the general problem is not solvable, we elaborate 
on the possibility to solve a restricted problem, searching 

5 Second kind Abel ODEs can be always rewritten as first kind 
Abel ODEs for w with the transformation go(x) + gi(x)y = 1/w. 
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for special solutions with H(—a) = H(—f3) = 0, i.e. with 
fi = in eq. (|38[) . This is done by writing, in full generality, 



H = k(T)(T + a) 2l (T + /3) 2m , 



(43) 



where / and m are integer numbers > 1, and k is a regular 
function at r = — a and r — —f5. The parity of the expo- 
nents forces k to be positive for r > —7, consistently with 
the non-negativity of H, and eq. (|38|) reduces to a linear 
ODE for /c 2 (r), so that only solutions everywhere positive 
are acceptable. Some algebra shows that the exact solution 



H± (r) = p(r) J k ± 2 / / (£) (t + 7) y/\t + a\\t + P\dt, (44) 



where fco is the free parameter. Note that the integral under 
square root is monotonic increasing (fo > 0), so that pos- 
it ivity is assured when ko > and the sign + is adopted. 
If the integral is convergent for r — )> 00, then also the sign 
— can be adopted, provided ko is larger than the value of 
the integral at infinity. In all cases, the square root cannot 
vanish for r — To > —7, as monotonicity of the integral then 
would produce negative values of k 2 for r > to. It follows 
that H diverges at r = —a and — f3 as p(r) (independently 
of the value of / and m), against the assumption, and re- 
duced solutions with fi do not exist. A simple argument 
shows that this negative it is to be expected. In fact, if one 
is allowed to fix H(-a) = H(-/3) = in eq. (|38|) . then the 
resulting non-homogeneous ODE becomes first-order linear 
for H 2 , and so it can solved in closed form. However, as 
the reduced equation is first order, its solution depends on 
a single parameter, and so in general only one of the two 
values H(—a) and H(—f3) can be set equal to zero. The ar- 
guments above show that the z-axis ODE for Vi potentials 
in dMOND is a genuine Abel equation of the second kind, 
to be solved numerically for assigned density profile. 

Concerning the asymptotic behavior of the density, fol- 
lowing the same treatment done for separable models in 
Sect. 3.1, it can be shown that also regular V\ models are 
characterized by a vanishing density at the center. The 
asymptotic analysis at infinity reveals some interesting dif- 
ference with the cases in Sect. 3.2. In fact, let us consider 
the Vi family 



F(t) = 



h(r) hiT 



h(r) 



1 for t 



(45) 



similar to the family of finite mass systems discussed in 
Sect. 3.2. We note that the parallel is only formal: while 
in the separable models the potential at large radii becomes 
spherical, in this case it remains ellipsoidal. This means that, 
at variance with the separable case, Vi models in the fam- 
ily above are characterized by infinite mass (or, in case of 
finite mass, by space sectors with negative density), because 



the potential of finite mass systems is dominated by the 
monopole term at large radii. 

We note that in Newtonian gravity eq. ([45]) with h = 1 
corresponds to the Binney (1981) triaxial logarithmic poten- 
tial, ln(l + x 2 /a 2 + y 2 /b 2 + z 2 /c 2 ) oc ln(A/xz/), fully discussed 
in ZP88 (Sects. 3.1 and 6.1 therein). Here we found that the 
Binney potential in dMOND leads to a non-spherical den- 
sity distribution of infinite mass, with a radial r~ 3 profile 
at large radii, with a negative density along the z-axis for 
sufficiently large t, independently of the values of a, b, c. In 
Cartesian coordinates: 

-2 + c 2 /a 2 + b 2 /a 2 



p(0,0,s) 



< c < b < a, 



00. (46) 



What happens if we allow for a more general h function in 
eq. (|45p ? For simplicity we do not repeat the analysis done 
in Sect. 3.2 for separable potentials, but we just report a few 
results. First, we found that the negative densities at large 
radii along the z-axis can be removed by appropriate choices 
of A and e in eq. ([27|) . but we were unable to construct 
everywhere positive mass models in the family (|45|) . Second, 
at variance with the separable case, at large radii p oc 1/r 3 
(modulated by an angular part) independently of e, so that 
finite mass dMOND systems in the V± family (|45|) do not 
exist. 

Obviously, other choices of F in eq. (|33|) are possible and 
worth of investigation. For example, ZP88 (eq. [2.13]) show 
that the Newtonian potential of the density distribution 

P * l + x 2 /a 2 + y 2 /b 2 + z 2 /c 2 ' (47) 

belongs to the Vi , and that the gradient of the potential can 
be written in terms of Carlson's (1979) symmetrized version 
of incomplete elliptic integrals. These functions are a closed 
family under differentiation, and therefore the dMOND den- 
sity distribution associated with this family can be written 
explicitely by using functions no more complicated than el- 
liptic integrals. Here, we do not study these models. 



6 SUMMARY AND CONCLUSIONS 

In this paper we studied the properties of density distri- 
butions obtained from MOND potentials separable in ellip- 
soidal coordinates. The investigation, besides the astrophys- 
ical context, is also interesting because the MOND operator, 
in the weak field limit (the so-called dMOND regime), re- 
duces to the non-linear p-Laplace operator V(|| V0|| p_2 V</>) 
with p — 3, the critical case in 3ft 3 . 

The main results obtained can be summarized as fol- 
lows: 

1. We proved that, as in the case of Newtonian grav- 
ity, also for MOND systems with separable potentials in 
ellipsoidal coordinates, the density profile along the z-axis 
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(the long axis of the ellipsoidal A-surfaces) determines the 
potential and so the density everywhere. The second-order 
ODE is however highly non- linear, and probably unsolvable 
in closed form, even in dMOND, in the limit of small flat- 
tenings and/or at large distances from the origin. Therefore, 
while we formally proved that the Kuzmin property holds in 
MOND, the associated Kuzmin formula is not known, and 
the Kuzmin Theorem remains unproven. 

2. However, we obtained rigorous asymptotic formulae 
for the density near the origin and at infinity, in case of reg- 
ular separable potentials. We showed that, at variance with 
the case of Newtonian separable systems, the density at the 
center vanishes, and we studied the order of vanishing as a 
function of the order of regularity of the potential. From this 
result it follows that MOND separable systems with regu- 
lar potentials are necessarily in the dMOND regime at their 
center, so that they would appear as centrally dark matter 
dominated (formally, with an infinite value for the mass- 
to-ligh ratio) when interpreted in the context of Netwonian 
gravity. Of course, this property may be shared by other 
sufficiently regular non-separable potentials in MOND. 

3. The analysis at large radii was performed under the 
assumptions of finite total mass system and a separable po- 
tential sufficiently regular to allow for a Frobenius expansion 
at infinity. In the regular case (i.e., when the expansion of 
the shape function h in eq. (|27|) reduces to a Taylor series 
in terms of integer powers of 1/r), the density at large radii 
is positive, provided a certain coefficient in the expansion is 
greater than a negative value (determined by the axial pa- 
rameters of the ellipsoidal coordinates). The density shape 
is not spherically symmetric, and the radial decline is pro- 
portional to ln(r)/r 5 , at variance with the Newtonian case 
of finite mass and finite flattening, when p oc 1/r 4 . In case 
of a genuine Frobenius expansion with exponent < e < 1, 
the total mass is still finite, but the density profile at large 
radii is instead spherically symmetric, with a milder radial 
decline p oc ln(r)/r 3+2e , and positive provided the constant 
mentioned above is positive. This resembles the Newtonian 
case of finite mass when the density profile along the sim- 
metry axis is steeper than 1/z 3 but shallower than 1/z 4 . 

4. We constructed some triaxial separable MOND mod- 
els, everywhere positive for all the explored axial ratios. The 
shapes range from almost perfectly ellipsoidal systems to 
curious systems with density depressions or overdensities 
(some of them similar to the models constructed in Ciotti 
et al. 2006), depending on the specific choice for the shape 
function h: these last models are unlikely to be useful for 
the description of real stellar systems. 

5. We briefly addressed the properties of the class of 
V\ potentials introduced in ZP88: these potentials, albeit 
non-separable, are simpler than the separable case, and they 



are known to obey the Kuzmin theorem in Newtonian grav- 
ity. We showed that they obey the Kuzmin property also in 
MOND, and we derived the ODE relating the potential to 
the density profile along the z-axis. This equation is consid- 
erably simpler than in the separable case, and in fact can 
be transformed in a Abel equation of second kind. Unfor- 
tunately, the general solution of this class of ODEs is not 
known, so the Kuzmin theorem cannot be proved. As an 
example of Vi systems, we studied the case of the MOND 
analogue of the Binney (1981) logaritmic potential, and we 
showed that the density profile becomes negative along the 
z-axis, while the radial profile declines as 1/r 3 . Some vari- 
ants of this potential however admits positive densities (at 
large radii) for some function /i, but still declining as 1/r 3 , 
and so being characterized by infinite total mass. 

6. Finally, we showed (Appendix B) that power-law ax- 
isymmetric potentials separable in parabolic coordinates, as- 
sociated with a central (weak) cusp can be constructed in 
MOND, in analogy with the family discovered by Sridhar & 
Touma (1997), albeit for a more restricted range of central 
slopes. If a central black hole is added, the central regions 
are still cuspy, but the cusp is Newtonian. 

We conclude by noting a few points that are rele- 
vant for successive investigations. The first is related to 
the phase-space distribution function. Presently our under- 
standing of the phase-space distribution function of MOND 
non-spherical systems still rely mostly on the numerical 
Schwarzschild method (Wang et al. 2008, Wu et al. 2009, 
2010) or N-body simulations (Nipoti et al. 2007ab, 2011), 
and the resulting systems are expected to be non-separable. 
Here we can derive some firm conclusion on MOND separa- 
ble models. In fact, the vanishing central density of triaxial 
regular models forces the associated distribution functions to 
vanish for the values of the three integral of motions allowing 
for orbits that cross the center. Now, the orbit classification 
for the Newtonian case assumes that the third derivative of 
F(t) is negative everywhere (e.g., Kuzmin 1973, Hunter & 
de Zeeuw 1992). This is the case for the reference model with 
F(t) — — r 2 ln(r)/2, leading to a third derivative which is 
in fact —1/r so indeed negative as we choose r > —7 > 0. 
Moreover, in case of a Frobenius h function, it is easy to 
prove that the expansion coefficients can be chosen so that 
negativity is assured (leaving true the positivity of the den- 
sity). In all these case, the models are supported by the 
four major orbit families. This indicates that these MOND 
models have the same orbit structure as Newtonian sepa- 
rable systems. If the condition is violated, then there are 
other/more orbit families possible. Restricting to the first 
cases, it makes plausible that selfconsistent models with van- 
ishing central density might well exist by populating the 
tube orbits only, leaving the box orbits out, as these would 
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be contributing positive density in the center. The machin- 
ery to construct such models with thin tubes only (zero ra- 
dial action) is available from Hunter & de Zeeuw (1992): by 
leaving the boxes out, this avoids having to compute the box 
orbit distribution function by numerically solving a large set 
of linear equations. The thin tube distribution functions are 
given in closed form, once F(r) is chosen and the density is 
known. 

The second point is that the obtained results would 
change if the fi function in eq. (3) has a non-zero lower 
bound, i.e. fi{t) —> /io for t — > 0. In fact, the deepest gravity 
regimes observationally probed in isolated galaxies are ~ 
O.Olao (e.g., see Zhao 2007, Famaey et al. 2007, Wu et al. 
2008). If such a lower limit for fi exists, then the dynamics 
would be Newtonian at very large radii and at very small 
radii, opening the possibility of MOND Stackel models with 
finite total mass and non-zero central density. 
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APPENDIX A: ELLIPSOIDAL COORDINATES 

We report here the main properties of ellipsoidal coordinates 
relevant for the present work; for a full discussion and further 
references on the subject, see Z85 and ZLB85. 

Al Definitions 

Ellipsoidal coordinates (A, /x, v) are curvilinear othogonal co- 
ordinates defined as the three real roots for r of the cubic 
equation 



+ 



y 



, = 1, (Al) 

r + a ' t + /3 r + 7 

where without loss of generality we assume a < f3 < 7 < 0, 
so that < —7 < v < —/3 < fi < —a < A. Surfaces of con- 
stant A are ellipsoids, surfaces of constant /i are hyperboloids 
of one sheet, and surfaces of constant v are hyperboloids of 
two sheetfl For very large values of A, the ellipsoidal sur- 
faces become more and more similar to spheres of radius 
\[X. The relations between (A, v) and the Cartesian coor- 
dinates (x,y, z) of a given point are 

2 _ ( A + a) (fi + a) (y + a) 



(a- P){a-i) 
2 (A + + + 
V (0_ a )(0_ 7 ) ' 

2 _ (A + 7 )(M + 7)(^ + 7) 



(A2) 



so that the origin corresponds to A = —a, \i — —f3 and v — 
—7. A series expansion shows that near the origin Cartesian 
and ellipsoidal coordinates are related by the (first order) 
asymptotic relations 



A + a, 



Z 2 rsj V -f 7. 

The metric coefficients are 
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6 As some confusion often arise on this point, we remark that the 
z axis is the long axis of the ellipsoidal A-surfaces, but usually it is 
the short axis of triaxial mass models with Stackel potential (e.g., 
see Fig. 1 in ZLB85). See also Section 4 in de Zeeuw et al. 1986). 



where 

a T = (r + a)(r + P)(t + 7); 



(A5) 



note that a a > 0, a v > 0, but a M < 0, consistent with the 
posit ivity of the metric coefficients. The gradient operator 
reads 



V — — — + ^-— - + — — - 

h\ dX hfj, dji hv dv ' 



(A6) 



where (ex^e^^e^) is the local basis of mutually orthogonal 
unitary vectors (e.g., Arfken & Weber 2005). Therefore, the 
squared norm of the gravitational field becomes 

2 



2 1 wa^ 2 

o x w, w v <9/i J hl\dv 

while the differential operator in eq. @ is 



1 dcj) d 



d 



+ 



h\ dX dX h\ d\i d\i hi dv dv ' 



v'. + v 2 + v 2 , 



Finally, the Laplace operator can be written as 
V 2 

where 

2 



(A-/x)(A-i/) 



d 2 da\ d 
'dA 5 + ~d\d\ 



(A7) 



(A8) 



(A9) 



(A10) 



and and V 2 , follow from the equation above by the rota- 



tion A — Y fi 



A, applied once and twice, respectively. 



A2 The leading terms of density expansion at the 
center 

With heavy but straightforward computation it can be 
shown that the coefficients Ai, Bi, and C\ appearing in 
eqs. (|19[) - (|20|) . needed in the density expansion near the cen- 
ter of generic separable dMOND system, are 

A 1 = £j^A 1 , Bi = ^Bi, C 1 = ^C 1 , (All) 

where A is defined in eq. (fl8l) . and 

( Ai = F(-a)(/3 + 7 - 2a)(/3 - 7) - F'(-a)A+ 
F(-/?)(a-7) 2 -F(- 7 )(a-/3) 2 , 
B 1 = F(-p)(a + 7 - 2/?)(a - 7) + F'{-/3)A+ 

F(-a)(/?- 7 ) 2 -F(-7)(«-/?) 2 , 1 ' 

d = F(- 7 )(a + p- 2 7 )(a - /?) - F'(- 7 )A+ 

F(-a)(/?- 7 ) 2 -F(-/?)(a- 7 ) 2 . 

In the special case Ai — B\ — Ci =0, the values of 
F' at the center are fixed by the vanishing of the system 
above. The coefficients of the resulting expansions reported 
in eq. (|22[) are 



A 2 



7 -A 2 , B 2 = ^^B 2 , C 2 = -^/c 2 ,(A13) 
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16 Ciotti, Zhao & de Zeeuw 



(A14) 



where now 
■ A 2 = F"(-a)A - 2F(-a)(P - 7)+ 

2F(-p)(a-j)-2F(- 1 )(a-p), 
B 2 = F"{-p)A - 2F(-a){/3 - 7)+ 

2F(-/3)(a - 7) - 2F(- 7 )(a - /?), 
C 2 = F"(- 7 )A - 2F(-a)(/3 - 7)+ 
2F(-/?)(a - 7) - 2F(- 7 )(a - f3). 
The additional request of regularity, i.e., A2 — B2 — 
C2 = 0, fixes the values of F" at the center from the van- 
ishing of the system above. The coefficients of the resulting 
expansions reported in eq. (|23|) are 

/3- 7 ; 



A 3 

B 3 ■- 
C 3 -- 



A 
a — 7 

A 

a-j3 
A~ 



-F"'(-a) 
F"'{-0), 
F'"(- 7 ). 



(A15) 



Finally, a similar analysis shows that the (y, z) plane is 
also separated in two regions, 

Sx = \{y,z): + -fl- < 1 1 , (A20) 

I p — a 7 — a J 

obtained by fixing A = —a in eq. (|A2[) , and the complemen- 
tary region = 3l 2 /S\, obtained for \i — —a. Solving the 
resulting systems and asking for positivity of the functions 
N = v + 7 and M = p + /3 (in 5a), or L = A + /3 (in 
we now get 



where 



(p(-a,M-P,N--y), (y,z)eSx, 
\p(L-p,-a,N->y), (y,z)eS», 



(A21) 



L — M — — , S -f-/3 > 0, (A22) 

R 2 = y 2 + z 2 , and N = R 2 -M. 



A3 Cartesian planes in ellipsoidal coordinates 

In the study of the shape and density distribution of triax- 
ial mass models expressed in ellipsoidal coordinates it may 
be helpful to have the expression for the coordinate planes 
0), (x,0, 2), and (0, y,z) in terms of the ellipsoidal 
coordinates. The following formulae can be easily deduced 
by simple geometrical arguments (see also de Zeeuw et al. 
1986). 

The (x,y) plane is obtained by requiring that z = in 
eq. (|A2[) . This request leads to fix v — —7, and solve for 
assigned x 2 and y 2 the resulting system for L = A + a > 
and M = fi + f3 > 0. Therefore, 



p(x,y,0) = p(L - a,M - P,->y), 



(A16) 



where the density at the l.h.s. is intended to be expressed in 
Cartesian coordinates, 



R 2 - S + ^{R 2 -S) 2 + A5x 2 



5>P-a>0, (A17) 



L, R 2 =x 2 + y 2 . 



and M — R 

The situation is slightly more complicated for the other 
two planes. In fact, the (x,z) plane is fully covered by the 
union of the region 



< 1 



(A18) 



1-P 13- 

obtained by fixing \i — —ft in eq. (|A2|) . with the complemen- 
tary region S v = 3ft 2 /S^, obtained for v — —ft. Solving the 
resulting systems and asking for positivity of the functions 
L = A + a and N = v + 7 (in S^), or M = \i + 7 (in S^), 
one gets 



p(x,0,2) 



-/3,iV- 7 ), (x,^)e^ 



\ p(L - a, M - 7, (x, z) G S u 

where L is again given by eq. ()A17|) but now R 2 
S = 7 - a > 0, and M = N = R 2 - L. 



(A19) 
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APPENDIX B: THE POWER-LAW 
AXISYMMETRIC SEPARABLE MODEL 

A property common to all triaxial systems with a Newtonian 
potential separable in ellipsoidal coordinates is a constant 
density core. In the axisymmetric case, Sridhar & Touma 
(1997) were able to construct separable potentials support- 
ing a central density cusp. The question is whether such 
property carries on in MOND as well. We show that this is 
in fact the case. We start from the separable potential 



(r+,r_) : 



[(l + cos6>) 3 - fe + (l-cos6>) 3 - fc ] 



(Bl) 



where the parabolic coordinates r + and r_ are related to 
the standard spherical coordinates by the identities 



r+ = r(l + cos 9), r_ = r(l — cos 6). 



(B2) 



Sridhar & Touma (1997) show that the density distribution 
associated with the potential above, via the Laplace opera- 
tor, is 

(2 - k - cos0)(l + cos6>) 2 - fc + (2 - k + cos0)(l - cos6>) 2 - fc 
P ^ , 

and discuss its properties as a function of k. In particular, for 
< k < 1 the density is cusped at the origin and positive 
everywhere, and so they conclude that cuspy systems (of 
infinite mass) with separable potential exist, at least in the 
axisymmetric case. In the critical case, k = 1, the density is 
cuspy, but it vanishes on the z-axis, being p oc sin 2 6/r. 

We do not embark on the interesting but long discussion 
of the density related to the potential (Bl) in MOND, but 
we note the following results. First, it is easy to show that 
the force depends on radius as r 1_fc , and does not vanish 
along any direction (as the radial component of V0 never 
vanishes). It follows that for k > 1 the MOND system will 
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Figure Al. Isodensity contours of the dMOND analogues of the Sridhar & Touma (1997) power-law axisymmetric cuspy models with 
separable potential in parabolic coordinates. From left to right, k = 1/3, 2/3, and 1. As explained in the test, the density diverges at the 
origin for 1/2 < k < 1, while it vanishes for < k < 1/2. For k = 1/2 the density (not shown) is independent of r, i.e., it is constant on 
cones with the common vertex on the origin. 



be similar - near the origin - to the Newtonian case, when 
however the density is unphysical (Sridhar & Touma 1997). 
For < k < 1 the system is in the dMOND regime near 
the origin (and Newtonian at infinity), while in the critical 
k = 1 case the force is independent of r, and so the system 
can be constructed in the dMOND regime everywhere. An 
application of the 3-Laplace operator to the family of power- 
law potential with < k < 1 easily shows that the density 
near the center (or everywhere, in the k — 1 case), can be 
written as 

pocr 1 - 2 */^), (B4) 

where the explicit function fk is easily calculated but is not 
reported here. Remarkably, in the range < k < 1 the func- 
tion fk is nowhere negative. From eq. (B4) it follows that the 
density vanishes at the center for < k < 1/2, is indepen- 
dent of r for k — 1/2 (but with different values along differ- 
ent directions), and it is cusped for 1/2 < k < 1. Moreover, 
fi oc sin 2 as in the Newtonian case. Therefore, separable 
potentials with a central (weak) cusp can be constructed also 
in MOND, but only in the restricted range 1/2 < k < 1. 
However, the resulting densities, albeit positive, are still 
characterized by an infinite total mass, and their shapes 
are quite unnatural, as can be seen from Fig. Al, where 
some examples are presented. Sridhar & Touma (1997) also 
showed that a black hole can be added at the center of these 
models, leaving separability unaffected. This remains true in 
MOND, as the gravitational field of the black hole switches 
the field near the center from dMOND to Newtonian, so that 
the system will be cuspy also for < k < 1/2: however, this 
is not a "genuine" MOND cusp. 

We finally note that the superposition of potentials of 
the family (Bl) with different values of k is still a separa- 
ble potential. Of course, the dMOND operator is non-linear, 



so that the associated densities are not the sum of the sep- 
arate components. However, we performed some numerical 
experiments, and we found that the resulting densities (for 
< k < 1) are still positive. This could open the way to 
the construction of new families of axisymmetric MOND 
systems with separable potentials and more general density 
distributions than pure power-laws. 



